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Abstract 

The off-lattice Boltzmann (OLB) method consists of numerical schemes which are used to solve the discrete 
Boltzmann equation. Unlike the commonly used lattice Boltzmann method, the spatial and time steps are 
uncoupled in the OLB method. In the currently proposed schemes, which can be broadly classified into 
Runge-Kutta-based and characteristics-based, the size of the time-step is limited due to numerical stability 
constraints. In this work, we systematically compare the numerical stability of the proposed schemes in terms 
of the maximum stable time-step. In line with the overall LB method, we investigate the available schemes 
where the advection approximation is explicit, and the collision approximation is either explicit or implicit. 
The comparison is done by implementing these schemes on benchmark incompressible flow problems such as 
Taylor vortex flow, Poiseuille flow and, lid-driven cavity flow. It is found that the characteristics-based OLB 
schemes are numerically more stable than the Runge-Kutta-based schemes. Additionally, we have observed 
that, with respect to time-step size, the scheme proposed by Bardow et al /? / is the most numerically 
stable and computationally efficient scheme compared to similar schemes, for the flow problems tested here. 
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1. Introduction 

The lattice Boltzmann (LB) method is an alternative and powerful numerical technique used for modeling 
a variety of complex hydrodynamic flows [? ? ]. Unlike conventional numerical methods which discretize 
the macroscale governing equations directly, the LB method solves a fully-discrete kinetic equation for 
distribution functions (DFs) fi(x, £), designed to reproduce the Navier-Stokes equation in the hydrodynamic 
limit. The LB method has advantages such as ease of parallelization, simplicity of programming, and a 
capability for incorporating model interactions for simulating complex flows. 

A defining feature of the LB method is the coupling between the velocity and space-time discretizations. 
That is, for a particular discrete-velocity set, £*, the coupling automatically fixes the temporal and spatial 
steps through the relation Ax = ^At. This procedure has some advantages such as numerical-diffusion free 
(exact) advection and computational efficiency (copy-operation). The coupling is, in fact, a carryover from 
the earliest LB models, which were based on Lattice Gas Automata (LG A). However, the LG A link was 
broken when it was shown more than a decade ago that the LB method can be derived directly from the 
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discrete Boltzmann equation as a special finite-difference scheme [???]. Consequently, the velocity-space 
can be discretized according to the flow-physics to be modeled. The discretization of space and time is a 
numerical requirement and, importantly, is not tied to the discretization of the velocity-space. 

As a consequence, a subset of the LB method, called the off-lattice Boltzmann (OLB) method, was 
developed where space and time are independently discretized, i.e. Ax ^ ^At. In the OLB method, we 
do not have the simplicity of a Lagrangian-type of evolution (streaming), rather the evolution of fa takes 
place in an Eulerian sense. The earliest OLB schemes were geared mainly towards extending the geometric 
flexibility of the LB method, which was previously limited, due to the requirement of a uniform Cartesian 
mesh. Several OLB schemes with different spatial discretization methods such as finite-volume (FV), finite- 
element (FE), and finite-difference (FD), along with their variants, have been developed. For example, OLB 
schemes were used for non-uniform mesh [? ], curvilinear co-ordinates [? ? ], unstructured mesh [???], 
finite element mesh [? ? ] among others. These advancements have made the LB method feasible for many 
practical engineering problems. 

In addition to improving the geometric flexibility of the LB method, OLB schemes can also be used to 
solve the discrete Boltzmann equation (DBE) with higher-order lattices. Higher-order lattices are sets of 
discrete velocities, which are more suited to model more complex flows such as thermal flows, micro-scale 
(high Knudsen number) flows, etc. In many of these velocity sets (also termed as non-space-filling or off- 
lattice), the discrete velocities cannot be expressed as an integer multiple of the smallest non-trivial speed. 
The D2Q16 velocity-set listed in [? ] and [? ] and the D2Q17 velocity-set in [? ] are typical examples. Since 
the regular stream-collide type of evolution scheme cannot be employed with these lattices, OLB schemes 
provide a viable evolution scheme for the DBE. 

While several sophisticated spatial-discretization methods have been developed, many of the studies use 
time-marching schemes such as explicit Euler or Runge-Kutta (RK) for temporal discretization. Typically, 
these schemes require very small values of At relative to the relaxation parameter r, to maintain numerical 
stability [? ? ]. This is in contrast to the LB method, which offers unconditional stability. Small At 
requirement is particularly restrictive in the case of flows with high Reynolds number (Re) flows where r 
is very small. Moreover, in the LB method, the Mach number Ma in the simulations has to be kept small 
(generally less than 0.1) to limit the compressibility errors. Small values of Ma lead to a slower convergence 
rate, especially for steady-state flow problems [? ? ]. Thus, the combined effects of small Ma and At increase 
the overall computational cost of the RK-based OLB schemes. 

Many alternative time-marching schemes have been proposed that maintain the numerical stability of 
the OLB method at higher values of At, relative to the relaxation parameter r, i.e. at higher At/r values. 
These schemes vary greatly in their numerical stability due to the different approximations of the collision 
and advection part of the DBE. Hence, there is a need to systematically compare their relative performance 
in terms of the numerical stability of these schemes. This work addresses this need. 

More specifically, we assess the stability of various OLB schemes, as quantified in terms of their maximum 
allowable At/r ratio. This is done via benchmark testing on incompressible flow problems such as Taylor- 
vortex flow, Poiseuille flow and lid-driven cavity flow. The on-lattice D2Q9 velocity set, which is used here 
for evaluation purposes, is described in Section 2.1. The various time-marching (OLB) schemes used in the 
comparative analysis are described in brief in Section 2.2. 
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2. Numerical Formulation 


2.1. Discrete Boltzmann Equation 

The basis for all OLB schemes is the Boltzmann equation with the Bhatnagar-Gross-Krook collision 
approximation [? ], which is given as: 

^+€-V/ = -l(/-/■*), (1) 

where / = /(#,£,£) is the single-particle distribution function, V/ = is the spatial gradient of /, £ 
is the microscale velocity, r is the relaxation time of the collision process, and f eq = / e<? (cc,£,£) is the 
local Maxwell-Boltzmann (equilibrium) distribution function. Equation 1 is continuous in velocity and 
configuration (x,t) space. To discretize the velocity space £, the equation is non-dimentionalized using a 
chosen speed of sound, and the resulting f eq is expanded in a Taylor-series of fluid velocity u up to second- 
order. The discrete velocities are then obtained from the requirement that the lower-order hydrodynamic 
moments with respect to the truncated f eq satisfy the conservation of mass, momentum, and energy [? ? ]. 
Following this procedure, we obtain the widely-used discrete velocity set of the D2Q9 lattice, for which the 
discrete Boltzmann-BGK equation can be written as: 

^+£i-v/ i = -t(/ i -/f), ( 2 ) 

where fi = /i(cc, t), f^ q = / z e<? (cc, t) and i = 0,1, 2 • • • ,8. Here, while the Greek subscripts a = {x, y} 
in 2D imply summation, the Latin subscripts (over velocity) do not imply summation. Equation 2 is termed 
as the discrete Boltzmann equation (DBE), and for a general class of discrete velocities, also referred to as 
the discrete velocity model (DVM). The D2Q9 velocity set is given by: 

{ (0,0) fori = 0 

V3c s cos(6^, sinOi) fori = 1,2,3,4 (3) 

V3c s (V^(cos0i, sin 6i)) fori = 5, 6, 7, 8 

where 0i = (i — l)7r/2 for i = 1 — 4, 0i = (2i — 9)7r/4 for i = 5 — 8, and c s = l/y/3 is the speed of sound in 
the lattice. Figure 1 shows a representation of the D2Q9 lattice. 



7 4 8 

Figure 1: D2Q9 velocity lattice. 


The discrete form of the equilibrium distribution function (EDF) is given by: 
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/f = PWi 
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( 4 ) 


i 9(£i • u ) 2 3m 2 
~ + 2c 3 2? 

for i = 0, 

fori = 1,2, 3,4, (5) 

fori = 5, 6, 7,8. 

The macroscale density and velocity are related to the DF through: 

8 8 

f = £/< = £/'" ( 6 ) 

i =0 i=0 

2 = 0 2 = 0 

It can be shown that a Chapman-Enskog expansion with the above discrete form of the EDF recovers the 
incompressible, isothermal Navier-Stokes equation in the limit of small Knudsen and Mach numbers with a 
shear viscosity v given by: 


+ 3 




where the weights, tOj, are: 


4/9 

Wi = \ 1/9 
1/36 


v = c 2 r (7) 

where r is the non-dimensional relaxation time. 

2.2. Off-Lattice Boltzmann Schemes 
2.2.1. Explicit Runge-Kutta based schemes 

Since the discrete velocities are constants, the DBE can be considered as a system of linear, first-order, 
ordinary differential equations (ODEs) with a weak source (collision) term. This assumption is generally 
valid only if the gradients of the conserved quantities in the flow are not too high, i.e., the collision term is not 
highly non-linear. The number of ODEs in the system equals the number of discrete velocities; for example; 
nine equations in case of the D2Q9 lattice. Therefore, in principle, commonly used time marching schemes 
for ODEs, such as Euler, Runge-Kutta, etc. can be employed for temporal discretization of Equation 1. On 
the other hand, FD or FV methods can be used for spatial discretization. 

Focusing on temporal discretization, a general second-order Runge-Kutta (RK2) based OLB scheme for 
Equation 2 can be written as: 


,n ^~ 2 _ 

2 

n At 

J i 2 

(8) 

pn+1 

2 

f? - 



where 

( 9 ) 

with the gradient term expanded in 2D Cartesian co-ordinates as: 
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df n 

£ • V f” = £ — 

^ Jz a dx a 



+ f,iy 



Here /™ +2 = fi(x,t n + ^), / 2 n+1 = fi{x,t n + At), etc. Similar expressions can be written for the 
fourth-order RK scheme (RK4) [? ]. In these schemes, the viscosity was related to the relaxation time 
through v m c 2 r. 

Many of the earliest OLB schemes employed the forward Euler, RK2 or RK4 schemes, in combination 
with a variety of spatial discretization schemes. Using the RK-based time marching schemes, the geomet¬ 
ric flexibility of the LB method was extended to non-Cartesian domains, non-uniform grids, body-fitted, 
stretched grids, etc. [?????]. In the case of FD spatial discretization, the spatial order-of-accuracy 
can also be increased arbitrarily, using higher-order Taylor approximations of the gradient terms. On the 
other hand, several discrete velocity models with non-space-filling velocity-sets also employed the RK-based 
schemes as the evolution equation [? ]. 

Despite the geometric flexibility made possible by the RK-based schemes, the size of the At relative to 
r has to be kept very small to maintain numerical stability. The constraint on At comes primarily from the 
explicit approximation of the advection (LHS) and collision (RHS) terms of the DBE. An explicit advection 
approximation imposes a stability criterion on the size of At through the CFL condition, CFL = At/Ax < 
1. The CFL condition is well-understood to be a necessary condition for the stability of advection type of 
equation. However, the overall stability of the scheme is governed by the more restrictive condition on At 
due to explicit approximation of collision, given by the approximate condition At < r [? ? ]. 


2.2.2. Characteristics-based schemes 

The characteristics-based OLB schemes are based on time integration of the DBE along the characteristics 
using the 0 —method [? ? ]: 


fn+l 

J i 


= j? + A t 


(i - 0)ft?) + o &} +1 


( 11 ) 


where we denote ff q ) for brevity; a tilde indicates a term on the characteristic line, i.e., 

fl 1+1 = fi(x + ^At,t + At), ft 1 = fi(x,t)] and 0 < 6 < 1. For 9 = {0, |, 1}, we obtain an explicit 0(At), 
an implicit 0(At 2 ), and an implicit O(At) approximation of the collision term, respectively. The O(At) 
schemes are Euler-type schemes, and the 0(At 2 ) is a Crank-Nicholson-type scheme. However, a 0(At 2 ) 
collision approximation can still be obtained for any value of 0 E [0,1], if the viscosity and relaxation time 
are related through: 


"=(xt~ 0 - 5+9 ) c2sAt - (12) 

A distinction should be made on the order of magnitude of r as used in the standard LBM versus in the OLB 
method. In the standard LBM, r is 0(1), and in fact typically in the range of 0.6 < r < 3.5. In the OLB 
method, on the other hand, due to non-dimentionalization, r is 0(Kn ), where Kn is the Knudsen number. 
Therefore, in OLB method, r is 0(0.01) or lower, depending upon the Re. 

Broadly based on Equation 11, several OLB schemes have been proposed that are numerical stable at 
much larger A t/r . In general, these schemes differ in their advection and collision approximations (explicit 
or implicit), as described below. 

Following Equation 11, Lee and Lin [? ] proposed a fully-explicit scheme (advection and collision), 
termed herein as AE/CE, which can be written as: 
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with A = v/(? s + 0.5 At corresponding to # = 0, where A is also termed as the modified-relaxation parameter. 

It is well-known in the ODE theory that an implicit approximation of the relaxation (collision) term is 
critical for numerical stability, especially for stiff equations (highly non-linear flows). Following this fact, 
a scheme similar to AE/CE but with an implicit collision (Cl) treatment i.e., 0 = 1/2 was proposed in [? 
], herein referred to as the AE/CI scheme. The AE/CI scheme, however, required an iterative predictor- 
corrector type of approximation for collision, which increases the computational cost and complexity of the 
simulations. 

In order to avoid an iterative procedure to approximate an implicit-collision, a variable transformation is 
often employed in the LB method to mask the implicitness of the collision term [? ? ]. A similar technique 
is employed in the OLB context, where a new DF, is defined as: 


9 i{x,t) = fi(x,t ) - A teCli(fi(x,t)). (14) 

Importantly, the variable transformation process preserves mass and momentum conservation, i.e., p = 
JA fi = gi and pu = JA ^ i f i = JA ^ i g i . The variable transformation with 0 = \ also maintains 0(At 2 ) 
accuracy of the collision approximation [? ]. 

With the variable transformation transformation technique and 0 = |, Guo and Zhao [? ] proposed an 
OLB scheme that can be written as: 

g n + l =f n_ M ,JJl _ _ f e q , ny (1 5) 

In this scheme, although the collision is implicit and 0(At 2 ), the advection is explicit and O(At) along the 
characteristics. Following Equation 12 for 0 = 1/2, the viscosity-relaxation time relation in this scheme is 
v = tc 2 s . In this work, this scheme is referred as GZ scheme. 

Bardow et al. [? ] combined the variable transformation technique, along with an explicit 0(At 2 ) 
advection approximation along the characteristics to yield the BKG scheme: 


n+1 _ 




c . +A _l !(„» _ Q eq ’ n ) 
^ a dx a + X (9t 9% ’ 


AC d 


r)n n 2 

t u yj I ( n n eg,n\ 


(16) 


At 3 d 

“2+ ia d+ 




d(gr-gf’ w ) 

dxij 


Due to the implicit treatment of the collision term through variable transformation, the At < r constraint no 
longer applies. The only constraint on At is the CFL condition due to explicit treatment of advection. The 
viscosity is related to the relaxation time through A = z//c 2 + 0.5At. A similar scheme called the unstructured 
lattice Boltzmann with memory (ULBEM) was proposed by [? ]. The various schemes described above are 
summarized in Table 1. 
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Table 1: Summary of various off-lattice Boltzmann schemes. Here E indicates explicit, and I indicates implicit. Spatial 
order-of-accuracy for advection depends upon the spatial discretization method and schemes used, and hence is excluded. 


Scheme 

Advection 

Collision 

v — A relation 

RK2/RK4 

0(At 2 ) E 

0{At 2 ) E 

A = v/6 2 s 

AE/CE 

OjAF) e 

ofAF) e 

X = i V c s + 0.5At 

AE/CI 

O(AC) E 

W( A?) i 


GZ 

O(At) E 

ojAF) i 

A = y~R s 

BKG 

0(At 1 ) E 

0(At*) I 

A = v jcf + 0.5At 


2.3. Implementation aspects of a OLB scheme 

Equation 16 represents a typical temporal evolution scheme for the DFs used in the OLB method. 
Clearly, the usual stream-collide scheme type of evolution on the standard LBM is replaced by a procedure 
that involves solving a set of ODEs at each node. For the BKG scheme, the algorithm consists of the following 
steps: 


1. Initialize ff and g\ according to prescribed initial conditions, i.e., set = g® = g* 9 ’ 0 , where g eq ^ = 
g^ q (p°,u°). p° and u° are the given initial conditions. 

2. For a particular time-step n, evaluate gf through Equation 14 (variable transformation). 

dq n d 2 n n 

3. Evaluate the gradient terms -^p-, etc. In the case of FD spatial discretization, with a central- 
differencing scheme, these quantities can be computed as: 


dg? _ gi{x + Ax,y,t n ) -gi(x- Ax,y,t n ) , A 2 > 
dx 2Ax + [ X h 

d 2 9i _ 9 i(x + Ax,y,t n )-2g i (x,y,t n )+g i (x-Ax,y,t n ) , 2 ^ 

dx 2 Ax 2 +U{/\x ). no 

Similar expressions can be written for the first-order derivatives in the ^/-direction, and the second-order 

a 2 

mixed derivatives -?[? ]. Other schemes such as upwind-differencing can also be employed. 

4. Evaluate per Equation 16. 

5. Evaluate p n+1 , u n+1 and g^ 9 ' n+1 . 

6. Evaluate / 2 n+1 per equation 14. 

7. If u n+1 has converged, then stop, if not, repeat steps 2-6. 

3. Numerical Tests 

To test the numerical stability of the various schemes, as a function of the maximum stable A t/r, several 
steady and unsteady, incompressible, two-dimensional flows were simulated, and their results are presented 
in this section. Since the focus of this work is on evaluation of the temporal schemes, the simpler finite- 
difference method is used to discretize the spatial domains. For uniformity, Equation 17 is used to evaluate 
the derivatives in the gradient term for all the OLB schemes tested here. The central-difference scheme is 
chosen since they are less diffusive than other 0( Ax 2 ) schemes. This minimizes the contribution of numerical 
diffusion to the overall stability of the scheme. The numerical stability of a scheme can be concluded by the 
maximum allowable At for a particular Re and grid size. In other words, for a particular flow problem, we 
fix the Re (fixed r) and grid size (fixed Ax) and vary At, until the simulation becomes unstable. 

The schemes were coded in C++ using the Armadillo linear algebra library [? ], and parallelized for a 
shared-memory architecture using OpenMP. 
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3.1. Taylor Vortex Flow 

To test the stability and accuracy of various schemes without the artifacts of a boundary treatment, we 
first simulate the Taylor-Green vortex flow. This flow represents the unsteady flow of a freely decaying two- 
dimensional vortex, and is often used to evaluate the effective viscosity and temporal and spatial accuracy of 
a scheme. Here, the flow is computed within a periodic square box defined as —7r < x, y < tt with a uniform 
Cartesian mesh of size 100 x 100 on the domain. The analytical solutions of the flow-field are given by: 


u{x,y,t) 
v(x,y,t ) 

p(x,y,t) 


u re f cos(kix) sin(fc 2 y) exp[— v{k\ + h%)t], 

k 2 2 )t], 


ki 

Urefy- sin(fcix) cos(k 2 y ) exp[— v{k\ 
K2 


Pref ^ 


k 2 

cos(2Aqic) + —| cos(2k 2 y) 

rCo 


exp[—2 v(k\ + k%)t\. 


(18) 


In our simulations, to minimize compressibility effects, the Mach number is set as 0.01. The reference 
velocity is u re f = Ma x c s , the reference density is p re f = 1, the reference pressure is p re f = p re //c^, and 
the wave numbers Aq and k 2 are constants set to 1.0 and 4.0, respectively. The non-equilibrium initialization 
scheme proposed by [? ] is used to initialize the DFs, and periodic boundary conditions are applied in the 
x and y directions. The Re of the flow is 100. 

Figures 2 and 3 show the horizontal and vertical velocity profiles at different times as obtained using 
the BKG scheme. The velocity profiles are shown for t/t c = 0.5,1, and 2, where t c = In 2 /v(k\ + k%) is the 
time at which the amplitude of the decay is halved. Importantly, these simulation have been obtained with 
At = 10r, where r = v/(? s . The corresponding CFL number is 0.7. This confirms the observations of high 
A t/r made in [? ]. Moreover, the average relative error as defined by: 


^ 5^7/ (I ^num U ex act \ T \Vnum Vexact \) 

E = —-- 

y (l^eccactl T |^eccact|) 

is < 1%. Here u nurn is the simulated velocity, u exact is the analytical velocity, and the summation is over 
vertical plane at x = 0. This result demonstrates that the BKG scheme successfully overcomes the At < r 
restriction imposed by the RK-based OLB schemes. The only remaining restriction on At is the local CFL 
number which is due to an explicit advection approximation. 
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Figure 2: Horizontal velocity profile of Taylor-Vortex flow 
at various times. The numerical results have been ob¬ 
tained with a At = lOr using the BKG scheme. 



Figure 3: Vertical velocity profile of Taylor-Vortex flow 
at various times with At = lOr using the BKG scheme. 


For comparison, we also simulated the Taylor-Green vortex flow using the GZ and the AE/CE schemes. 
The GZ scheme, although having treated the collision term implicitly, could not achieve higher A t/r values. 
In fact, in their simulations, A t/r ~ 1.6 for a Re ~ 63 with a mixed (central-upwind) differencing scheme. 
Furthermore, for a purely central-difference scheme, i.e., without the additional marginal stability of upwind 
schemes, the authors showed analytically using the von-Neumann stability analysis, that the maximum A t/r 
that can be attained is less than 1.5, and that at low CFL numbers. We have observed that the AE/CE 
permits At > r, but stable simulations are obtained at much lower values of CFL number (CFL < 0.2) 
compared to the BKG AE/CI scheme. 

3.2. 2-D Plane Poiseuille Flow 

Another useful numerical test is Poiseuille flow. A plane Poiseuille flow describes the steady, laminar flow 
of an incompressible fluid in a rectangular channel driven by a pressure gradient. Assuming symmetry and 
incompressibility, it can be shown that for Poiseuille flow, the Navier-Stokes momentum equation reduces 
to: 


d 2 u dp 
^dy 2 dx 

which has a exact steady state solution for velocity given by: 


(19) 


u(y) = 4u max j^(l-^ for 0 <y<H (20) 

v(x,y) = 0 

where, H is the channel height and is the center-line velocity where the magnitude of velocity is the 

maximum. The Reynolds number of the flow is defined by Re = u max H/u. In a LB simulation, imposing a 
pressure difference by specifying inlet and outlet densities increases the compressibility errors [? ]. Therefore, 
the effect of A p is imposed on the flow through an equivalent body-force F. This force has the same effect as 
having a pressure-gradient in the channel which produces the chosen u max . This body force can be evaluated 
from F = 8u max pu/H 2 . 
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For the reference case, the BKG scheme is used to obtain the results presented here. In the simulation, 
we set H = 1 and L = 1, where L is the channel length. The domain is discretized into a uniform Cartesian 
mesh of size 100 x 100. The Mach number based on u max is set to 0.1732, which corresponds to a u max of 
0.1. The velocity is initialized to zero everywhere and the average density is set to one. Since a steady-state 
flow is simulated, the equilibrium initial condition is applied to the DFs. A non-equilibrium extrapolation 
boundary scheme is applied on the no-slip top and bottom walls [? ], and a periodic boundary condition is 
applied at the inlet and outlets. Various Re flows are simulated by varying the viscosity but keeping the Ma 
constant. The simulations are run until a steady state criterion, defined as: 


C = 




- U, 


n—50 | 

i,j I 


Ei 


< 10 " 


\ u i,j\ 


is attained. Here u 7 ^ - = u(xi,yj,nAt), where n is the time-step number, and the summation is over the 
entire flow field. 

The simulated steady-state stream-wise velocity profile for Re = 100, along with the analytical solution 
is shown in Figure 4. Clearly, the simulated and the analytical values are in excellent agreement with each 
other. The inset shows the velocity profile close to the wall, which is also in close agreement to analytical 
values. This serves as validation for the applicability and accuracy of the non-equilibrium boundary condition 
used in the simulation. It is worth noting that the velocity profile was generated on a 100 x 100 grid with a 
At = 5r. Even at such a high At, the average relative error as defined earlier is ^ 2%. 



y/L 


Figure 4: Steady-state horizontal velocity profile for Poiseuille flow at Re = 100. The results were obtained using At = 5 r on 
a 100 X 100 mesh using the BKG scheme. 


For comparison, the plane Poiseuille flow is also simulated using the AE/CE scheme, and the GZ scheme. 


10 










Figure 5 shows the maximum allowable CFL number versus Re for the steady Poiseuille flow for the different 
schemes. This plot is obtained by fixing the grid size to 100 x 100 for each scheme, and then varying At. 
This is repeated for the various Re, as shown in the figure. 

It can be seen that at low Re , both the AE/CE and BKG schemes have almost similar maximum CFL 
numbers. However, as Re increases, the maximum allowable CFL number is clearly much higher with 
the BKG scheme, as compared to the AE/CE scheme. This enhanced stability at higher Re results from 
the implicit treatment of the collision term, which becomes highly non-linear as Re increases, thus directly 
benefiting from the local implicit treatment. Importantly, the implicit treatment of collision, allows for 
At > r, thus resulting in large At and thereby reducing the computational effort. 



Figure 5: Maximum allowable CFL number vs. Re for steady Poiseuille flow. 


3.3. Lid-Driven Cavity Flow 

Isothermal, 2-D lid-driven cavity flow is commonly used as a benchmark case to test and evaluate nu¬ 
merical schemes for incompressible viscous flows. The flow domain consists of a square cavity with three 
stationary walls on the sides and bottom, and a top wall (lid) that moves with a uniform tangential ve¬ 
locity uud. Despite the simple geometry, this flow exhibits many complex flow patterns such as formation 
of vortices near corners due to singularities. In our context, the lid-driven cavity case is also useful for 
quantitatively evaluating the effects of collision approximations (explicit/implicit) on flows with moderately 
large non-linearities (high Re ), and thus quantifying the overall stability of different schemes. 

The computational domain consists of a square with height L = H = 1. The top wall moves with a 
constant velocity of uud = 0.1, and the Reynolds number of the flow is Re = uudH/u. The domain is 
discretized on a uniform mesh of size 257 x 257. The non-equilibrium extrapolation scheme is applied for 
all of the walls. The flow is initialized by setting p = 1 and u = 0 in the entire flow. The flow is simulated 
using four schemes: BKG, AE/CE, GZ, and also a RK2-based scheme. 

As before, for the reference case, we present the results of simulations using the BKG scheme. Steady-state 
velocity components along the horizontal and vertical centerlines at various Re are shown in Figures 6 and 7. 
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These results have been obtained with At = 4 r and CFL = 0.5. The results are compared with results from 
Ghia et al. who obtained their results with a 257 x 257 grid using the coupled strongly implicit multi-grid 
method, and a vorticity-stream-function formulation [? ] . The velocity profiles change from curved at 
lower Re to linear for higher Re , which is consistent with the benchmark solutions. The near-linear velocity 
profiles at higher Re in the central core of the cavity indicates a region of uniform vorticity. 



Figure 6: Centerline horizontal-velocity profile for a lid- 
driven cavity at various Re. The solid lines indicate the 
simulated values, while the markers (■, A, ▼) indicate the 
corresponding solution from Ghia et al. [? ] 



Figure 7: Centerline vertical-velocity profile for a lid- 
driven cavity at various Re. The solid lines indicate the 
simulated values, while the markers (■, A, ▼) indicate the 
corresponding solution from Ghia et al. [? ] 


To assess and quantify the effects of the size of At on numerical stability, in Figures 8-13, we plot the 
stability region for each scheme, as determined by two parameters: At/r and At/Ax (please note that the 
scales of the axes vary across the figures). Broadly speaking, while At/r represents the non-linear stability 
constraint imposed by collision, At/Ax represent the linear-stability constraints of explicit advection. At/Ax 
also represents the CFL number in the case of a D2Q9 lattice. Hence, a map determined by these two 
parameters should serve a guide to gauge the overall numerical stability of the schemes with respect to size 
of At. In all of the maps, • indicates a stable solution and o indicates an unstable solution. Figure 8 shows 
the stability region for the RK2-based OLB scheme with central differencing for Re = 100. 
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Figure 8: Stability region for the RK2 scheme at Re = 

100. Figure 9: Stability region for the GZ scheme at Re = 100. 

From Figure 8, we can observe that the RK2-based scheme is unstable beyond A t/r > 2, irrespective 
of the CFL number. Moreover, the stable solutions which are obtained at A t/r < 2 are done so at very 
low CFL values. As described before, the small At requirement stems from the fact that as Re increases, 
the collision term becomes more non-linear and stiff, which requires an implicit approximation for numerical 
stability. Stability of RK2-based schemes can be increased slightly by adopting second-order upwind schemes; 
however, the marginal stability is added at the cost of increasing numerical diffusion. Similarly, multi-stage 
Runge-Kutta schemes such as RK4 can also increase stability. However, all RK-based OLB schemes involve 
multiple evaluation of the f^ q term for each advancement in At, the number of evaluations depending upon 
the stages in the scheme. Since evaluation of f eq is a computationally intensive step, RK-based schemes 
are also computationally inefficient. Therefore, we can conclude that RK2 based OLB schemes are not 
particularly suitable for flows with large non-linearities. 

Figure 9 shows the stability region for the GZ scheme for Re = 100. From the stability-region plot, 
we can again see that the scheme is unstable for all A t/r > 2, irrespective of the CFL number. Thus, in 
spite of the implicit collision approximation, large At could not be used in the scheme. This observation is 
consistent with the stability analysis presented by the authors of the scheme. 

The stability regions of the AE/CE scheme are shown in Figure 10 and 11. From the figures, it is evident 
that the AE/CE scheme does allow A t/r > 2, but does so only at lower CFL number (CFL < 0.2). As 
Re increases, however, the explicit collision approximation becomes inadequate, and hence the At has to be 
smaller. 
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Figure 10: Stability region for the AE/CE scheme at Re = n c , ..... . , ,. A _, , D 

7 Figure 11: Stability region for the AE/CE scheme at Re = 

Finally, Figures 12 and 13 show the stability regions for the BKG scheme at Re = 400 and 1000. We 
can observe that the simulations are stable at A t/r as high as 30, and that at high CFL numbers. The 
trend also does not deteriorate with increasing Re in the range that we have tested. This demonstrates the 
unconditional collision stability of the BKG scheme, with At restricted only by the local CFL number due 
to explicit advection. 

The BKG scheme is also computationally more efficient than the comparable AE/CI scheme. This is 
because whereas a predictor-corrector type of scheme is needed for the implicit collision approximation in 
the AE/CI scheme, a simple variable transformation given by Equation 14 is required in the BKG scheme. 
Additionally, the BKG scheme, as with all advection-explicit OLB schemes, retains the data-locality feature 
of the LB method, and hence can be easily parallelized. 

Finally, it is important to note that we have discussed only the numerical stability of different OLB 
schemes in terms of maximum allowable At. However, as for all explicit advection schemes, including the 
BKG scheme, much smaller values of Ax and At are required to get accurate solutions. 
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Figure 12: Stability region for the BKG scheme at Re = 
400. 


Figure 13: Stability region for the BKG scheme at Re 
1000. 
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4. Conclusions 


In this work, several explicit OLB schemes have been compared by implementing them for benchmark 
flow problems. The following conclusions can be drawn from the observations: 

1. The characteristics-based OLB schemes provide higher numerical stability compared to RK-based 
schemes. 

2. In characteristics based-schemes, the At < r constraint no longer applies, even at high Re. 

3. The scheme proposed by Bardow et al. is stable over a much wider range of simulation parameters 
(At/r, At/Ax) compared to other similar characteristics-based OLB schemes for the problems tested 
in this work. The BKG scheme also retains the simple explicit form of the LB method, while providing 
unconditional collision-stability. 

These conclusions indicate that the BKG scheme provides the most stable and efficient explicit time-marching 
scheme, which can be extended to flow problems with FV or FD discretization. This scheme, in theory, can 
also be adopted for thermal problems with both off- and on-lattice discrete-velocity sets [? ]. 
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